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Abstract 

We discovered a class of self-similar solutions in nonlinear models describing the for- 
mation of morphogen gradients, the concentration fields of molecules acting as spatial 
regulators of cell differention in developing tissues. These models account for diffu- 
sion and self-induced degration of locally produced chemical signals. When production 
starts, the signal concentration is equal to zero throughout the system. We found that 
in the limit of infinitely large signal production strength the solution of this problem is 
given by the product of the steady state concentration profile and a function of the dif- 
fusion similarity variable. We derived a nonlinear boundary value problem satisfied by 
this function and used a variational approach to prove that this problem has a unique 
solution in a natural setting. Using the asymptotic behavior of the solutions established 
by the analysis, we constructed these solutions numerically by the shooting method. 
Finally, we demonstrated that the obtained solutions may be easily approximated by 
simple analytical expressions, thus providing an accurate global characterization of the 
dynamics in an important class of non-linear models of morphogen gradient formation. 
Our results illustrate the power of analytical approaches to studying nonlinear models 
of biophysical processes. 



1 Introduction 

Reaction-diffusion processes are involved in multiple aspects of embryogenesis. In partic- 
ular, a combination of extracellular diffusion and degradation of locally produced proteins 
can establish concentration fields of chemical signals that control spatial and temporal 
gene expression patterns in developing tissues [ij. Such concentration fields are known as 
morphogen gradients and have been identified in contexts as diverse as neural development 
in vertebrates and wing morphogenesis in insects (2l[3]. 
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A canonical model of morphogen gradient formation is given by the following initial 
boundary value problem pHsl: 



8C rflC 

= ^ 9^ - C{x,t = Q) = Q, (1) 

dC 



dx 



= Q, C{x = oo,t) = 0. (2) 

x=0 



Here C = C{x, t) is the concentration of a morphogen as a function of distance x > to the 
tissue boundary and time t > 0. The morphogen is produced with a constant rate Q at the 
tissue boundary (x = 0), diffuses with diffusivity D in the tissue (x > 0) and is degraded 
in the tissue following some rate law characterized by the pseudo first-order rate constant 
k(C) > 0. This model provides a minimal description of complex biochemical and cellular 
processes in real tissues, and has been recently used to quantitatively describe morphogen 



gradients in a number of experimental systems [4|[9 -11 . 

The level of gene expression in a particular cell located a certain distance away from the 
signal source depends only on the local concentration of that signal (or, more generally, on 
its time history at that location) . Some of the genes controlled by morphogen gradients are 
directly involved in regulating processes of cell differentiation. Another subset of genes act 
more indirectly, regulating morphogens themselves. For example, morphogens can control 
the expression of cell surface molecules that accelerate the rate of morphogen degradation, 
establishing a feedback loop. Such self-induced morphogen degradation can be modeled by 
having the degradation constant k{C) to be an increasing function of concentration. 

The presence of cooperative effects that are commonly observed in the transcriptional 
responses to morphogen gradients [S] suggests to consider a general class of degradation 
rates given by a power law |4|: 

k{C) = knC''-\ n>l. (3) 

This type of nonlinearity was first considered by Eldar et al. , who demonstrated that power 
law degradation kinetics may generate morphogen gradients that are robust with respect 
to large variations in the source strength [2]. They based their conclusions on the analysis 
of the steady version of Eq. ([T|. Specifically, they demonstrated that, unlike the solutions 
of the corresponding linear problem, i.e., Eqs. ([T]) and ^ with k{C) = const (in which 
the solution depends on Q multiplicatively) , the stationary solution Cs{x) of Eqs. ([l])-([3]) 
approaches an asymptotic limit when Q — )• oo. As a consequence, the steady state of a 
system operating in the regime of large Q's will be insensitive to variations in the strength of 
the source. This has important implications for robustness of steady morphogen gradients 
established by localized production, diffusion and self-induced degradation. 

In this paper we show that robustness of the steady state solutions of Eqs. Q-Q dis- 
cussed above carries over to the solutions of the full time-dependent problem. Remarkably, 
we found that for large values of Q the solution of the initial boundary value problem given 
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by Eqs. ([T])-([3]) approaches a self-similar form: 



C{x,t) = Cs{x)^{x/VDt), (4) 

where is a universal function of ^ = x/\/T)t which depends only on n and decreases 
monotonically from i;^)=lat.^ = 0tO(/> = 0at^ = oo. The self-similar profile function 
(f){(^) is obtained by considering the singular version of the initial boundary value problem 
with Q = oo. 

Our results can be summarized as follows. We derived a nonlinear boundary value 
problem satisfied by (/){£,)■ We used a variational approach to prove that this problem has 
a unique solution in a natural setting. Using the asymptotic behavior of the solutions 
established by the analysis, we constructed these solutions numerically by the shooting 
method. We also showed that the obtained solutions may be easily fitted to simple ana- 
lytical expressions, thus providing an accurate global characterization of the dynamics in 
an important class of non-linear models of morphogen gradient formation. For example, in 
the biophysically important case n = 2, in which the feedback is mediated by the simplest 
bimolecular interaction our approach yields 

» , "t^^I^ ^. (5) 

4000 + ^9 
4000 + 5^663^ 



- -77: 1,2 ' ™ = 2- (6) 



A comparison of the prediction of Eqs. ^ and Q with the numerical solution of Eqs. 
([T])-([3]) for a representative set of parameters is presented in Fig. [l| 

This paper is organized as follows. We first introduce the non-dimensionalization of 
the considered initial boundary value problem and give scaling arguments that lead us 
to consider singular self-similar solutions of the problem. We then undertake a rigorous 
analysis of existence, uniqueness and qualitative properties of the considered self-similar 
solutions, which are summarized in Theorem 1. After that, we perform a numerical study 
of self-similar solutions whose existence was established in the Analysis section and discuss 
their implications to the dynamics. Finally, we discuss the obtained results in connection 
with the studies of other types of self-similar solutions in the considered problem, their 
significance for the threshold crossing and conclude with some open problems. 

2 Results 

In this section, we formulate a general setting for the analysis of Eqs. ([T])-([3]) and present 
our main findings. 
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Figure 1: Comparison of the numerical solution of Eqs. ([T])-|3]) for n = 2,L' = l, A;2 = l 
and Q = 10^ (thin solid lines) with the prediction of Eqs. ([5| and ([g]) (dashed lines) at 
times t = 0.04, 0.16, 0.64, 2.56, 10.24. Thick solid line shows the steady state. 



2.1 Non-dimensionalization 

We begin by introducing the following dimensionless variables 



x' = |, t' = |, n=^, (7) 



where 



L=i t] , T= ^, (8) 



and Co is some reference morphogen concentration, corresponding, e.g., to the threshold 
of expression of a downstream regulated gene. In these new variables, the initial boundary 
value problem in Eqs. Q-Q takes the form 

Ut = Uxx - (x, t) G [0, oo) X (0, oo), 

n^(0,t) = -a tG(0,oo), (9) 
u(x,0) = X G [0,oo). 



where 



a = ^ ^ (10) 



is the dimensionless source strength. 
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2.2 Scaling arguments 

Let us now discuss the approach of the solutions of Eq. ^ to the unique steady state, 
which for this problem is given explicitly by the following expression [12j : 

. ^ f 2(n+l) 1 ^ 

Vaix) = . (11) 

[ [(2"(n+ l)ai-")"+i + (n - l)x] J 

It is not difficult to see that u{x, t) — )■ Va{x) from below as t — )• oo, implying that the fraction 
of the steady concentration u(x, t)/va{x) reached at a given point x > at time i > will 
approach unity for t ^ 1 [12]. In view of the diffusive nature of the processes involved 
in establishing the steady concentration profile, one may expect that the approach to the 
steady state occurs on the scale associated with diffusion. Therefore, to better understand 
the dynamics, we plot this fraction versus x/y/i for the solution of Eq. ^ with n = 2 and 
a = 1 obtained numerically for several values of t. The result is presented in Fig. [2} One 
can see from Fig. [2] that the solution of Eq. ^ at different values of t collapses onto a 
single master curve for t ^ 1. Furthermore, increasing the value of a makes this collapse 
sooner. We also checked that the same phenomenon occurs for different values of n. This 
strongly suggests [l3] the existence of a hidden self-similarity in the underlying dynamical 
behavior of the solutions of Eq. ^ . 

We note that the solutions of Eq. ^ are invariant with respect to the following scaling 
transformation: 

a = Xa, t = X i+" t, X =X"+^x, u = X"+'^u. (12) 
In other words, increasing the source strength a by a factor of A decreases the time scale 

2(n-l) 

of approach to the steady state by a factor of A "+i at fixed value of xjyt. Therefore, 
the approach to the universal curve in Fig. [2] must occur on the time scale 

2(1-") , , 

~ a "+i . (13) 

This scale was recently identified by us in the analysis of the local accumulation time in 
the particular case of Eq. ([9]) [12]. Observe that ^ as a — )• oo for all n > 1. Thus, 



our numerical results suggest that in the limit a — )• oo the ratio u{x,t)/va{x) depends only 
on xl\ft for all t > 0, exhibiting self-similar behavior. 



2.3 Singular solutions and the similarity ansatz 

The numerical observations discussed in the preceding section suggest the need to consider 
the following singular initial boundary value problem: 

ut = Uxx-u^ (a;,t) G (0,oo) X (0,oo), 

n(0,t) = oo, te(0,oo), (14) 
u{x,0) = X £ (0,oo). 
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Figure 2: Example of the collapse of the solutions of Eq. ^ onto a universal master curve 
at large times. Results of the numerical solution of Eq. B with n = 2 and a = 1. Thin 
lines show snapshots of the solution corresponding to t = 0.1, 1, 10, 100 (the direction of 
time increase is indicated by the arrow). The bold line shows the asymptotic master curve. 



By a solution to Eq. (14), we mean a classical solution for all {x,t) E (0,oo) x (0,oo) 
decaying sufficiently fast as x — )• +oo for all t > 0, and continuous up to t = for all x > 0. 
Note that for each n > 1 this problem possesses a singular stationary solution 



2(n + l) 



1 2 
n— 1 / \ n — l 



(15) 



(n — 1)2 J \x ^ 

which is the limit of Va{x) as a — t- oo for each x > 0. Therefore, in view of the discussion 
above, solution of Eq. ( 14 ) is expected to take the form 

u{x,t) = Voo{x)(l){x/Vi), (16) 

for some universal function with values between zero and one, which depends only on 
n. 



We now substitute the similarity ansatz from Eq. (16) into Eq. (14). After some 



algebra, this leads to the following equation for the self-similar profile 



4C 



n 



2(n + l) 
1 J ' (n- 1)2' 



+ 



[1 



0, 



(17) 



which must hold for all ^ G (0,oo). Consistent with the interpretation of Eq. (14), this 
equation needs to be supplemented with the boundary-like conditions 



lim 0(^) = 1, lim (/)(^) = 0. 



(18) 
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Existence and multiplicity of solutions of Eqs. ( 17 ) , ( 18 ) are not a priori obvious in view 



of both the non-linearity and the presence of singular terms in the considered boundary 
value problem. Below we establish existence and uniqueness of these solutions for all n > 1 
within a natural class of functions. Later, in the following section, we construct these 
solutions numerically. 



2.4 Analysis 

This section is concerned with the proof of existence and uniqueness of solutions of Eqs. 



(17), n8| in a natural setting. The reader interested in the application of our analysis 



to Eqs. may skip this part and proceed directly to the following section, which 

discusses the numerical solution of Eqs. (17), (18). 



c 



For the purposes of the analysis it is convenient to rewrite Eq. (17), using a new variable 
= In^: 



+ 



,2C 



n + 3 



n 



dcj) 2(n + 1) 



1)2 



:i 



m-l^ 



0, 



where (j) S C'^iJR.) and has the respective limits 



lim m 



lim 0(C) 



(19) 



(20) 



We will prove existence and uniqueness of solutions of Eqs. ( 19 ) and ( 20 ) in the weighted 



Sobolev space -ff^(M, d/u), which is obtained as the completion of the family of smooth 
functions with compact support with respect to the Sobolev norm ||-||j7i(M,d^), defined as 



m(R,dii) 



Ikcll 



L2(R,dfi) 



+ \w\ 



|2 

Il2( 



f^w'^{()dij,{(), and the measure dfi is 

2C /n + 3 



p{C)dC, p{C) = exp 



C 



(21) 



(22) 



4 \n- 1 

Our existence and uniqueness result is given by the following theorem. 

Theorem 1. There exists a unique weak solution (p of Eq. (19), such that (p — r] € 
H^{]S.,diJ.), with fi defined in Eq. (22), for every rj £ C°°(M), such that r]{C) = 1 for 
all C < and r/(C) = for all ( > 1. Furthermore, (p G C°°(M), satisfies Eq. ^ 



classically and < (p < 1. In addition, (j) is strictly decreasing and satisfies Eq. ( 20 ) . 
Proof. The proof consists of five steps. 



Step 1. We first note that Eq. ( 19 ) is the Euler-Lagrange equation for the energy functional 




+ 



T] 



n 



1 



'{n + 1- 2(^' 
(n-l)2 



>dfj,, 



(23) 
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where r]{Q is as in the statement of the theorem. The natural admissible class A for £ is: 



1] G H^{R,dfi), < < 1}. 



(24) 



Note that the role of r] in the definition of £ is to ensure that the integral in Eq. (23) 
converges for all (j) £ A. The precise form of rj{Q is unimportant. 

Step 2. We now establish weak sequential lower-semicontinuity and coercivity of the func- 
tional £ in the admissible class A in the following sense: let (/>fc = r/ + w^, where Wk ^ w 
in H^{M.,dfi). Then 1) lim ini /^^^ £((pi:) > £{(!)), where (j) = rj + w, and 2) if £{(t)k) < M 
for some M G M, then | < M' for some M' > 0. 

Let us introduce the notation <?(</>, {a^h)) for the integral in Eq. (23), in which Integra- 

(25) 



tion is over all C, G (a, b). Arguing as in [14[ Lemma 4.1], for ^ 1 we have 



R 



wldiJL yw £ H^{R,dfi). 



Therefore, from Eq. (25) we find that 

£{^k,{R,+oc)) > fe2^ 



n + 1 



wld^, 



(26) 



V (n-l)2 

which is positive for i? S> 1. Similarly, taking into account that the integrand in Eq. 
(23) is non-negative for C < 0, we have £{(j)k, {—oo,0)) > 0. Since <S(-,(0,i?)) is lower- 
semicontinuous by standard theory 15 , we obtain £{(j)k) > £{4'k, (0, R)), yielding the first 
claim in view of arbitrariness of i? S> 1. 



To prove coercivity, we first note that by Eq. ( 25 ) 
£{cl)k,{R, +oo)) > [ 

JR 



1 / dwk 

2 V ^^C 



(n 



"+1 2I , 



> 



1 



dwk\' 
dC J 



+ wl } d^i, 



(27) 



for i? » 1. On the other hand, since n — I — cjp'in + 1 — 2(/)" ^) > (n — 1)(1 — 4')'^ for ah 
< < 1, we have 



f(</>fc,(-oo,0)) > 



1 / dwk\ 
2\dC ) 



+ 



n 



dfj,. 



(28) 



Finally, by boundedness of (j)k and r], we also have 



£i<PkAo,R)) > 



R 



dwk 
dC 



+ wl > dfi — CR, 



(29) 
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for some C > 0, whenever | |_ffi(R,d/i) is sufficiently large. So the second claim follows. 
Step 3. In view of the lower-semicontinuity and coercivity of E proved in Step 2, by the 



direct method of calculus of variations there exists a minimizer 
(t!) = and 6 



1 solve Eq. (19), we also have that 



^ AoiE. Noting that 
is a weak solution of Eq. (19) by 



continuous differentiability oiE in H^{E.,dfi). Furthermore, by standard elliptic regularity 
[16 



theory 16 



and is, in fact, a classical solution of Eq 
, we have < cp < 1. To show monotonicity, suppose 



(19). Also, by strong 
to the 



maximum principle 

contrary, that (/>(a) < <f){b) for some a < b. Then 4>{C) attains a local minimum for some 
Co £ (—00,6). However, by Eq. (19) we have d? 4>{C,q) / dC,"^ < 0, giving a contradiction. By 
the same argument d(p/dC, = is also impossible for any C G M. Finally, since (p — rj € 



H^(K,dfi), monotonicity implies Eq. (20). 

Step 4. We now discuss the asymptotic behavior of the minimizers obtained in Step 3 as 
C — )■ +00. Performing the Liouville transformation by introducing ip = G L'^{R^ +00)5 
where p is defined in Eq. (22) and i? G M is arbitrary, we rewrite Eq. (19) in the form 



Here q{C) = qo{C) + qi{C), where 

qoiC) = 
QiiC) = 



C>R- 



1 

4 

2(n + 1 



e^^ n ■ 
4 n 



=2C 



+ 1 



(30) 

(31) 
(32) 



Observe that q{C) > qo{C) > | > for all C > ^ > 1- Therefore, Eq. ^ has two 
linearly- independent positive solutions ipi and ip2, such that tpi ^ and "02 — co together 

^ CV'i G L2(i?,+oo) for 



with their derivatives as ( 
some C > 0, and 



+00 (see e., 



17 ). In particular, 



1 — rt 

qiiO = o{p 2 ), 



(33) 



so qiiC) has a super-exponential decay as C — +00. 

Now let ipQ be the unique positive solution of Eq. ( pO| ) with q = qo and ipoiR) = 1 
which goes to zero as C — )• +00. Then we claim that V'l (0/^^0(0 ~^ c for some < c < cxd. 
Indeed, after straightforward algebraic manipulations we have 



d 
dC 



In 



qi{s) , , , , ds. 



c 



(34) 



The result then follows by the decay of ipQ and ipi and the estimate in Eq. (33) upon 



integration of Eq. ( 34 ) . 
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Step 5. We now prove uniqueness of the obtained solution, taking advantage of a sort of 
convexity of E pointed out in |18|. Suppose, to the contrary, that there are two functions 
(t>i,4>2 ^ -A which solve Eq. (19) weakly. Define = ^/t(j)^~+~(Y^~^^^ ^ since in view 



of the result of Step 4 we have m < (f)i / (f)2 < M for some M > m > 0. It is easy to see that 
the function E{t) := ^(0*) is twice continuously differentiable for all t G [0, 1]. A direct 
computation yields 




blf{t<t>l + {l-t)4>l)^ U/i(C). (35) 



Therefore, d'^E{t)/dt^ > for all t G [0,1], and so E{t) is strictly convex. However, 
since the map t i— )•(/>* — r/ is of class C"'^([0, 1]; /7^(M, d/x)), this contradicts the fact that 



dE{0)/dt = dE{l)/dt = by the assumption that 4'i and (j)2 solve weakly Eq. (19) and 
hence are critical points oi 8. □ 

Remark 1. Note that the arguments in Step 4 above imply that the asymptotic behavior 
of the self-similar profile (p{S,) as ^ ^ oo is the same as that of the decaying solution of Eq. 



(17) linearized around (f) = 0. A similar argument shows that the behavior of (/){£,) as ^ — )• 



is the same as that of the corresponding solution of Eq. (17) linearized around (j) = 1. 



2.5 Numerics 

We now construct the self-similar profiles, whose existence and uniqueness was established 
in Theorem 1, numerically for several values of n > 1. We use the shooting method to 



construct the solutions of Eq. (17), which requires knowledge of the asymptotic behavior 



of 0(^) near ^ = and ^ = oo. To obtain this behavior, we linearize Eq. (17) around 
the equilibria = and (j) = 1, which is justified by Remark [TJ Denote the corresponding 
solutions of the linearized equations as (pQ and (pi , respectively. By a direct computation 



^o(6 = Ci^^M 



1 1 e 



n-1'2' 



.C.e-n^.(l.^.i.f). (30) 

where M(a, b, z) and U{a, b, z) are the confluent hypergeometric functions of the first and 
second kind, respectively jl9|. Using the asymptotic expansions of these functions for large 



19 , one can see that (/>o(C) ~^ &s ^ — )• oo, if and only if the constant Ci = 0. Therefore, 



from the asymptotic expansion of U we have 



1 ^2 5 — n 

HO ~ e"4« , oo. (37) 
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Figure 3: Self-similar profiles for different values of n. Results of the numerical 
solution of Eqs. (17) and ( |18| ) for n = 1.25, 1.5, 2, 3, 4, 6, oo. The thick line is the graph of 
the function given by Eq. |m overlaying the profile for n = 2. 



Similarly 



2(n + l) 

MO = ? "-^ 



CiM 



+C2U 




(38) 



Once again, for a bounded solution at ^ = we must set C2 = 0, which leads to 

2(n+l) 

l-</'(0~C^^> C^O. (39) 

The results of the numerical solution of Eq. \17\ whose asymptotic behavior is governed 
by Eqs. (36) and (38) are presented in Fig. ^ One can see that the self-similar profiles 



form a monotonically decreasing family of functions parametrized by n. The solutions 
approach = 1 on finite intervals as n — >• 1 and </>(^) = 1 — erf(^/2) as n — ?• 00 (the 
latter solves Eq. (17) corresponding to n = 00). We also found that for the biophysically 
important case n = 2 the self-similar profile can be approximated by the simple expression 
given by Eq. ^ within ~ 1% accuracy. The graph of this function, which essentially 



coincides with that of the numerical solution of Eq. (17) is shown in Fig. [3] with a thick 
line. Note that this profile also coincides with the limiting profile in Fig. [2] for t 



00. 



2.6 Dynamics 



We now discuss the dynamical behavior of the obtained self-similar solutions of Eq. (14). 



The picture remains qualitatively the same for all n > 1, so in the following we restrict our 
attention to the biophysically important case of n = 2. 
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Figure 4: Self-similar solutions u{x,t) of Eq. (14) at several values of x for n = 2 



First consider the time course of the solution u{x,t) given by Eq. (16) at a fixed 



location, i.e. at a fixed value of x > 0. From the self-similarity ansatz in Eq. (16) it 
is clear that the time scale of these dynamics is governed by diffusion, i.e., t ~ A 
convenient characterization of local dynamical time scale can be made in terms of the 
local accumulation time r°°(x) = tp{x, t)dt, where the probability density-like quantity 

du(x,t) 



p{x,t) 



integration by parts, one oBIains 



12 



20 



Upon substitution of Eq. ( 16 ) into this formula and an 



r-(x) 



ax 



(40) 



where numerically a ~ 0.122. We note that by Eq. (39) the integral in Eq. (40) converges 



for all n > 1. The solution for several values of x is shown in Fig. l4| Furthermore, as 

^2 

follows from Eqs. (37) and (39), when t <^ t°^{x), we have u{x,t) ~ (x/t^'^)e~ 47 ^ which is 
exponentially small. At the same time, for t ^ t^(x) we have (^00(2;) — u{x,t)) /voo{x) ~ 
{T°°{x)/t)^, i.e., u approaches the stationary solution, with the distance to the stationary 
solution decaying as 0{t~^). 



We now consider the motion of the level sets of the solutions of Eq. (14). For a given 
c > 0, let us define Xc{t) as the unique value of x, such that u{x, t) = c for each t > 0. As 



follows from Eqs. (16), the function Xcif) can be determined parametrically as 



xc = mi)/c)^'\ t = 6m/ice), n 



(41) 



The graphs of Xc{t) for a few values of c are shown in Fig. [sj Once again, the dynamics of 
Xc can be characterized by the local accumulation time (xjf ) given by Eq. ( 40 ) , where 
: (6/c)^/^ is the asymptotic value of Xc{t) as t — )• cxd. One can see from Eqs. (37) and 



(41) that for t ^ r°°(x^) we have Xc — 2{tlnt ^)^/^. Thus, all level sets move together 
for short times, as can also be seen from Fig. [5j On the other hand, for t ^> t°°(x^) the 
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Figure 5: The positions Xc{t) of level sets {u{x,t) = c} of the self-similar solution of Eq. 



(14) at several values of c for n = 2. 



level set position Xc{t) approaches [t] = 0{t Within ~ 2% accuracy the 

functions Xc{t) can be approximated by the following simple expression: 

Mt)^f^"f-?5'°"'V''. « = 2, (42) 



1 + 0.76ct J 

This formula implies that Xc{t) comes within 5% of x^ at t ~ 2r°°(x^). 

3 Discussion 

We have provided an analytical characterization of the dynamics of morphogen concentra- 
tion profiles in models with self-induced morphogen degradation. Our results reveal the 
presence of self-similarity in the course of the approach of the concentration profiles to 
their steady states in either the limit of large source strengths or for large distances away 
from the source. In addition to demonstrating the self-similar nature of the dynamics, 
we rigorously established existence and uniqueness of the associated singular self-similar 
solutions to Eq. ([T]) and constructed these solutions numerically for several values of n. 
The obtained solutions may be readily used to study various characteristics of the local 



kinetics of morphogen concentration. In particular, Eqs. (41) and (42) obtained from 
the numerical self-similar solutions provide a characterization of threshold crossing events, 
which determine the times at which a morphogen gradient switches the gene expression on 
or off at a given point. 

3.1 Comparison with self-similar transients 

Mathematically, we have constructed a new class of self-similar solutions to Eqs. ([T]), ([s]) 
on half-line. These solutions can be trivially extended to the whole real line by a reflection 
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and can be viewed as singular solutions of Eqs. ([T]), ([s]) that blow up at the origin. We note 
that Eqs. ([T]), Q and their d-dimensional analog is known to possess another kind of self 
similar solutions which were extensively studied l21l|23|, starting from the works of |24,25 



(see also |26| for a less technical introduction and 27 for a variational treatment). These 
are classical solutions of Eqs. ([T]), ([s]) for t > 0, which concentrate to a point mass at t = 
and describe the transient dynamics for such initial data. They also play an important 
role in the long-time behavior of the associated initial value problem (see e. g. 24,28-31]). 
The solutions found by us can be viewed as the counterparts of the very singular solutions 
constructed in [25]. Our solutions are more singular than those of [25] in the sense that 
the singularity in the former is concentrated on a half-line (x = 0, t > 0) in the (x, t) plane, 
while the singularity in the latter occurs only at a single point (x = 0,t = 0). 



3.2 Extensions and open problems 

It would be interesting to understand the role our self-similar solutions play for the singular 
solutions of the initial value problem associated withEqs. ([T]) , ([sj) for general non-zero initial 
data. Let us point out that even the basic questions of existence and uniqueness of such 
singular solutions for the considered parabolic problems in suitable function classes are 



currently open (see 32 for a very recent related work). Other natural extensions include 
higher dimensional versions of the considered problem, as well as a proof of global stability 
of self-similar solutions. From the point of view of applications, it is also important to 
consider singular solutions of Eq. ([ij) for time-varying sources, for which both types of 
self-similar solutions that are present in the system may be relevant. 
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